home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 11 / Cream of the Crop 11-1.iso / math / ast51src.zip / CALC.C < prev    next >
C/C++ Source or Header  |  1995-12-31  |  23KB  |  717 lines
  1. /*
  2. ** Astrolog (Version 5.10) File: calc.c
  3. **
  4. ** IMPORTANT NOTICE: The graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1995 by Walter D. Pullen
  6. ** (Astara@msn.com). Permission is granted to freely use and
  7. ** distribute these routines provided one doesn't sell, restrict, or
  8. ** profit from them in any way. Modification is allowed provided these
  9. ** notices remain with any altered or edited versions of the program.
  10. **
  11. ** The main planetary calculation routines used in this program have
  12. ** been Copyrighted and the core of this program is basically a
  13. ** conversion to C of the routines created by James Neely as listed in
  14. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  15. ** available from Matrix Software. The copyright gives us permission to
  16. ** use the routines for personal use but not to sell them or profit from
  17. ** them in any way.
  18. **
  19. ** The PostScript code within the core graphics routines are programmed
  20. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  21. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  22. **
  23. ** The extended accurate ephemeris databases and formulas are from the
  24. ** calculation routines in the program "Placalc" and are programmed and
  25. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  26. ** (alois@azur.ch). The use of that source code is subject to
  27. ** regulations made by Astrodienst Zurich, and the code is not in the
  28. ** public domain. This copyright notice must not be changed or removed
  29. ** by any user of this program.
  30. **
  31. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  32. ** X Window graphics initially programmed 10/23-29/1991.
  33. ** PostScript graphics initially programmed 11/29-30/1992.
  34. ** Last code change made 12/27/1995.
  35. */
  36.  
  37. #include "astrolog.h"
  38.  
  39.  
  40. /*
  41. ******************************************************************************
  42. ** House Cusp Calculations.
  43. ******************************************************************************
  44. */
  45.  
  46. /* This is a subprocedure of ComputeInHouses(). Given a zodiac position,  */
  47. /* return which of the twelve houses it falls in. Remember that a special */
  48. /* check has to be done for the house that spans 0 degrees Aries.         */
  49.  
  50. int HousePlaceIn(rDeg)
  51. real rDeg;
  52. {
  53.   int i = 0;
  54.  
  55.   rDeg = Mod(rDeg + 0.5/60.0/60.0);
  56.   do {
  57.     i++;
  58.   } while (!(i >= cSign ||
  59.       (rDeg >= house[i] && rDeg < house[Mod12(i+1)]) ||
  60.       (house[i] > house[Mod12(i+1)] &&
  61.       (rDeg >= house[i] || rDeg < house[Mod12(i+1)]))));
  62.   return i;
  63. }
  64.  
  65.  
  66. /* For each object in the chart, determine what house it belongs in. */
  67.  
  68. void ComputeInHouses()
  69. {
  70.   int i;
  71.  
  72.   for (i = 1; i <= cObj; i++)
  73.     inhouse[i] = HousePlaceIn(planet[i]);
  74. }
  75.  
  76.  
  77. /* This house system is just like the Equal system except that we start */
  78. /* our 12 equal segments from the Midheaven instead of the Ascendant.   */
  79.  
  80. void HouseEqualMidheaven()
  81. {
  82.   int i;
  83.  
  84.   for (i = 1; i <= cSign; i++)
  85.     house[i] = Mod(MC-270.0+30.0*(real)(i-1));
  86. }
  87.  
  88.  
  89. /* Compute the cusp positions using the Alcabitius house system. */
  90.  
  91. void HouseAlcabitius()
  92. {
  93.   real rDecl, rSda, rSna, r;
  94.   int i;
  95.  
  96.   rDecl = RAsin(RSin(OB) * RSinD(Asc));
  97.   r = -RTan(AA) * RTan(rDecl);
  98.   rSda = DFromR(RAcos(r));
  99.   rSna = rDegHalf - rSda;
  100.   house[sLib] = DFromR(RA) - rSna;
  101.   house[sSco] = DFromR(RA) - rSna*2.0/3.0;
  102.   house[sSag] = DFromR(RA) - rSna/3.0;
  103.   house[sCap] = DFromR(RA);
  104.   house[sAqu] = DFromR(RA) + rSda/3.0;
  105.   house[sPis] = DFromR(RA) + rSda*2.0/3.0;
  106.   for (i = sLib; i <= sPis; i++)
  107.     house[i] = Mod(house[i]+is.rSid);
  108.   for (i = sAri; i <= sVir; i++)
  109.     house[i] = Mod(house[i+6]+rDegHalf);
  110. }
  111.  
  112.  
  113. /* This is a new house system similar in philosophy to Porphyry houses.   */
  114. /* Instead of just trisecting the difference in each quadrant, we do a    */
  115. /* smooth sinusoidal distribution of the difference around all the cusps. */
  116.  
  117. void HousePorphyryNeo()
  118. {
  119.   real delta;
  120.   int i;
  121.  
  122.   delta = (MinDistance(MC, Asc) - rDegQuad)/4.0;
  123.   house[sLib] = Mod(Asc+rDegHalf); house[sCap] = MC;
  124.   house[sAqu] = Mod(house[sCap] + 30.0 + delta   + is.rSid);
  125.   house[sPis] = Mod(house[sAqu] + 30.0 + delta*2 + is.rSid);
  126.   house[sSag] = Mod(house[sCap] - 30.0 + delta   + is.rSid);
  127.   house[sSco] = Mod(house[sSag] - 30.0 + delta*2 + is.rSid);
  128.   for (i = sAri; i < sLib; i++)
  129.     house[i] = Mod(house[i+6]-rDegHalf);
  130. }
  131.  
  132.  
  133. /* The "Whole" house system is like the Equal system with 30 degree houses, */
  134. /* where the 1st house starts at zero degrees of the sign of the Ascendant. */
  135.  
  136. void HouseWhole()
  137. {
  138.   int i;
  139.  
  140.   for (i = 1; i <= cSign; i++)
  141.     house[i] = Mod((SFromZ(Asc)-1)*30+ZFromS(i));
  142. }
  143.  
  144.  
  145. /* In "null" houses, the cusps are always fixed to start at their cor-    */
  146. /* responding sign, i.e. the 1st house is always at 0 degrees Aries, etc. */
  147.  
  148. void HouseNull()
  149. {
  150.   int i;
  151.  
  152.   for (i = 1; i <= cSign; i++)
  153.     house[i] = Mod(ZFromS(i));
  154. }
  155.  
  156.  
  157. /* Calculate the house cusp positions, using the specified algorithm. */
  158.  
  159. void ComputeHouses(housesystem)
  160. int housesystem;
  161. {
  162.   char sz[cchSzDef];
  163.  
  164.   if (RAbs(AA) > RFromD(rDegQuad-rAxis) && housesystem < 2) {
  165.     sprintf(sz,
  166.       "The %s system of houses is not defined at extreme latitudes.",
  167.       szSystem[housesystem]);
  168.     PrintError(sz);
  169.     Terminate(tcFatal);
  170.   }
  171.   switch (housesystem) {
  172.   case  1: HouseKoch();           break;
  173.   case  2: HouseEqual();          break;
  174.   case  3: HouseCampanus();       break;
  175.   case  4: HouseMeridian();       break;
  176.   case  5: HouseRegiomontanus();  break;
  177.   case  6: HousePorphyry();       break;
  178.   case  7: HouseMorinus();        break;
  179.   case  8: HouseTopocentric();    break;
  180.   case  9: HouseAlcabitius();     break;
  181.   case 10: HouseEqualMidheaven(); break;
  182.   case 11: HousePorphyryNeo();    break;
  183.   case 12: HouseWhole();          break;
  184.   case 13: HouseNull();           break;
  185.   default: HousePlacidus();
  186.   }
  187. }
  188.  
  189.  
  190. /*
  191. ******************************************************************************
  192. ** Star Position Calculations.
  193. ******************************************************************************
  194. */
  195.  
  196. /* This is used by the chart calculation routine to calculate the positions */
  197. /* of the fixed stars. Since the stars don't move in the sky over time,     */
  198. /* getting their positions is mostly just reading info from an array and    */
  199. /* converting it to the correct reference frame. However, we have to add    */
  200. /* in the correct precession for the tropical zodiac, and sort the final    */
  201. /* index list based on what order the stars are supposed to be printed in.  */
  202.  
  203. void ComputeStars(SD)
  204. real SD;
  205. {
  206.   int i, j;
  207.   real x, y, z;
  208.  
  209.   /* Read in star positions. */
  210.  
  211.   for (i = 1; i <= cStar; i++) {
  212.     x = rStarData[i*6-6]; y = rStarData[i*6-5]; z = rStarData[i*6-4];
  213.     planet[oNorm+i] = RFromD(x*rDegMax/24.0+y*15.0/60.0+z*0.25/60.0);
  214.     x = rStarData[i*6-3]; y = rStarData[i*6-2]; z = rStarData[i*6-1];
  215.     planetalt[oNorm+i] = RFromD(x+y/60.0+z/60.0/60.0);
  216.     /* Convert to ecliptic zodiac coordinates. */
  217.     EquToEcl(&planet[oNorm+i], &planetalt[oNorm+i]);
  218.     planet[oNorm+i] = Mod(DFromR(planet[oNorm+i])+rEpoch2000+SD);
  219.     planetalt[oNorm+i] = DFromR(planetalt[oNorm+i]);
  220.     ret[oNorm+i] = RFromD(rDegMax/26000.0/365.25);
  221.     starname[i] = i;
  222.   }
  223.  
  224.   /* Sort the index list if -Uz, -Ul, -Un, or -Ub switch in effect. */
  225.  
  226.   if (us.nStar > 1) for (i = 2; i <= cStar; i++) {
  227.     j = i-1;
  228.  
  229.     /* Compare star names for -Un switch. */
  230.  
  231.     if (us.nStar == 'n') while (j > 0 && NCompareSz(
  232.       szObjName[oNorm+starname[j]], szObjName[oNorm+starname[j+1]]) > 0) {
  233.       SwapN(starname[j], starname[j+1]);
  234.       j--;
  235.  
  236.     /* Compare star brightnesses for -Ub switch. */
  237.  
  238.     } else if (us.nStar == 'b') while (j > 0 &&
  239.       rStarBright[starname[j]] > rStarBright[starname[j+1]]) {
  240.       SwapN(starname[j], starname[j+1]);
  241.       j--;
  242.  
  243.     /* Compare star zodiac locations for -Uz switch. */
  244.  
  245.     } else if (us.nStar == 'z') while (j > 0 &&
  246.       planet[oNorm+starname[j]] > planet[oNorm+starname[j+1]]) {
  247.       SwapN(starname[j], starname[j+1]);
  248.       j--;
  249.  
  250.     /* Compare star declinations for -Ul switch. */
  251.  
  252.     } else if (us.nStar == 'l') while (j > 0 &&
  253.       planetalt[oNorm+starname[j]] < planetalt[oNorm+starname[j+1]]) {
  254.       SwapN(starname[j], starname[j+1]);
  255.       j--;
  256.     }
  257.   }
  258. }
  259.  
  260.  
  261. /*
  262. ******************************************************************************
  263. ** Chart Calculation.
  264. ******************************************************************************
  265. */
  266.  
  267. /* Given a zodiac degree, transform it into its Decan sign, where each    */
  268. /* sign is trisected into the three signs of its element. For example,    */
  269. /* 1 Aries -> 3 Aries, 10 Leo -> 0 Sagittarius, 25 Sagittarius -> 15 Leo. */
  270.  
  271. real Decan(deg)
  272. real deg;
  273. {
  274.   int sign;
  275.   real unit;
  276.  
  277.   sign = SFromZ(deg);
  278.   unit = deg - ZFromS(sign);
  279.   sign = Mod12(sign + 4*((int)RFloor(unit/10.0)));
  280.   unit = (unit - RFloor(unit/10.0)*10.0)*3.0;
  281.   return ZFromS(sign)+unit;
  282. }
  283.  
  284.  
  285. /* Transform spherical to rectangular coordinates in x, y, z. */
  286.  
  287. void SphToRec(r, azi, alt, rx, ry, rz)
  288. real r, azi, alt, *rx, *ry, *rz;
  289. {
  290.   real rT;
  291.  
  292.   *rz = r *RSinD(alt);
  293.   rT  = r *RCosD(alt);
  294.   *rx = rT*RCosD(azi);
  295.   *ry = rT*RSinD(azi);
  296. }
  297.  
  298.  
  299. #ifdef PLACALC
  300. /* Compute the positions of the planets at a certain time using the Placalc */
  301. /* accurate formulas and ephemeris. This will superseed the Matrix routine  */
  302. /* values and is only called with the -b switch is in effect. Not all       */
  303. /* objects or modes are available using this, but some additional values    */
  304. /* such as Moon and Node velocities not available without -b are. (This is  */
  305. /* the one place in Astrolog which calls the Placalc package functions.)    */
  306.  
  307. void ComputePlacalc(t)
  308. real t;
  309. {
  310.   int i;
  311.   real r1, r2, r3, r4;
  312.  
  313.   /* We can compute the positions of Sun through Pluto, Chiron, the four */
  314.   /* asteroids, Lilith, and the (true or mean) North Node using Placalc. */
  315.   /* The other objects must be done elsewhere.                           */
  316.  
  317.   for (i = oSun; i <= oLil; i++) {
  318.     if ((ignore[i] && i > oMoo) ||
  319.       (us.fPlacalcAst && FBetween(i, oCer, oVes)))
  320.       continue;
  321.     if (FPlacalcPlanet(i, t*36525.0+2415020.0, us.objCenter != oSun,
  322.       &r1, &r2, &r3, &r4)) {
  323.  
  324.       /* Note that this can't compute charts with central planets other */
  325.       /* than the Sun or Earth or relative velocities in current state. */
  326.  
  327.       planet[i]    = Mod(r1 + is.rSid);
  328.       planetalt[i] = r2;
  329.       ret[i]       = RFromD(r3);
  330.  
  331.       /* Compute x,y,z coordinates from azimuth, altitude, and distance. */
  332.  
  333.       SphToRec(r4, planet[i], planetalt[i],
  334.         &spacex[i], &spacey[i], &spacez[i]);
  335.     }
  336.   }
  337. }
  338. #endif
  339.  
  340.  
  341. /* This is probably the main routine in all of Astrolog. It generates a   */
  342. /* chart, calculating the positions of all the celestial bodies and house */
  343. /* cusps, based on the current chart information, and saves them for use  */
  344. /* by any of the display routines.                                        */
  345.  
  346. real CastChart(fDate)
  347. bool fDate;
  348. {
  349.   CI ci;
  350.   real housetemp[cSign+1], Off = 0.0, vtx, j;
  351.   int i, k;
  352.  
  353.   /* Hack: Time zone +/-24 means to have the time of day be in Local Mean */
  354.   /* Time (LMT). This is done by making the time zone value reflect the   */
  355.   /* logical offset from GMT as indicated by the chart's longitude value. */
  356.  
  357.   if (RAbs(ZZ) == 24.0)
  358.     ZZ = DegToDec(DecToDeg(OO)/15.0);
  359.   ci = ciCore;
  360.  
  361.   if (MM == -1) {
  362.  
  363.     /* Hack: If month is negative, then we know chart was read in through a  */
  364.     /* -o0 position file, so the planet positions are already in the arrays. */
  365.  
  366.     MC = planet[oMC]; Asc = planet[oAsc];
  367.   } else {
  368.     for (i = 1; i <= cObj; i++) {
  369.       planet[i] = planetalt[i] = 0.0;    /* On ecliptic unless we say so.  */
  370.       ret[i] = 1.0;                      /* Direct until we say otherwise. */
  371.     }
  372.     Off = ProcessInput(fDate);
  373.     ComputeVariables(&vtx);
  374.     if (us.fGeodetic)               /* Check for -G geodetic chart. */
  375.       RA = RFromD(Mod(-OO));
  376.     MC  = CuspMidheaven();          /* Calculate our Ascendant & Midheaven. */
  377.     Asc = CuspAscendant();
  378.     ComputeHouses(us.nHouseSystem); /* Go calculate house cusps. */
  379.  
  380.     /* Go calculate planet, Moon, and North Node positions. */
  381.  
  382.     ComputePlanets();
  383.     if (!ignore[oMoo] || !ignore[oNod] || !ignore[oSou] || !ignore[oFor]) {
  384.       ComputeLunar(&planet[oMoo], &planetalt[oMoo],
  385.         &planet[oNod], &planetalt[oNod]);
  386.       ret[oNod] = -1.0;
  387.     }
  388.  
  389.     /* Compute more accurate ephemeris positions for certain objects. */
  390.  
  391. #ifdef PLACALC
  392.     if (us.fPlacalc)
  393.       ComputePlacalc(T);
  394. #endif
  395.     if (!us.fPlacalc) {
  396.       planet[oSou] = Mod(planet[oNod]+rDegHalf);
  397.       ret[oSou] = ret[oNod] = RFromD(-0.053);
  398.       ret[oMoo] = RFromD(12.5);
  399.     }
  400.  
  401.     /* Calculate position of Part of Fortune. */
  402.  
  403.     j = planet[oMoo]-planet[oSun];
  404.     if (us.nArabicNight < 0)
  405.       neg(j);
  406.     j = RAbs(j) < rDegQuad ? j : j - RSgn(j)*rDegMax;
  407.     planet[oFor] = Mod(j+Asc);
  408.  
  409.     /* Fill in "planet" positions corresponding to house cusps. */
  410.  
  411.     planet[oVtx] = vtx; planet[oEP] = CuspEastPoint();
  412.     for (i = 1; i <= cSign; i++)
  413.       planet[cuspLo + i - 1] = house[i];
  414.     if (!us.fHouseAngle) {
  415.       planet[oAsc] = Asc; planet[oMC] = MC;
  416.       planet[oDes] = Mod(Asc + rDegHalf); planet[oNad] = Mod(MC + rDegHalf);
  417.     }
  418.     for (i = oFor; i <= cuspHi; i++)
  419.       ret[i] = RFromD(rDegMax);
  420.   }
  421.  
  422.   /* Go calculate star positions if -U switch in effect. */
  423.  
  424.   if (us.nStar)
  425.     ComputeStars(us.fSidereal ? 0.0 : -Off);
  426.  
  427.   /* Transform ecliptic to equatorial coordinates if -sr in effect. */
  428.  
  429.   if (us.fEquator)
  430.     for (i = 1; i <= cObj; i++) if (!ignore[i]) {
  431.       planet[i]    = RFromD(Tropical(planet[i]));
  432.       planetalt[i] = RFromD(planetalt[i]);
  433.       EclToEqu(&planet[i], &planetalt[i]);
  434.       planet[i]    = DFromR(planet[i]);
  435.       planetalt[i] = DFromR(planetalt[i]);
  436.     }
  437.  
  438.   /* Now, we may have to modify the base positions we calculated above */
  439.   /* based on what type of chart we are generating.                    */
  440.  
  441.   if (us.fProgress && us.fSolarArc) {  /* Are we doing -p0 solar arc chart? */
  442.     for (i = 1; i <= cObj; i++)
  443.       planet[i] = Mod(planet[i] + (is.JDp - is.JD) / us.rProgDay);
  444.     for (i = 1; i <= cSign; i++)
  445.       house[i]  = Mod(house[i]  + (is.JDp - is.JD) / us.rProgDay);
  446.     }
  447.   if (us.nHarmonic > 1)            /* Are we doing a -x harmonic chart?     */
  448.     for (i = 1; i <= cObj; i++)
  449.       planet[i] = Mod(planet[i] * (real)us.nHarmonic);
  450.   if (us.objOnAsc) {
  451.     if (us.objOnAsc > 0)           /* Is -1 put on Ascendant in effect?     */
  452.       j = planet[us.objOnAsc]-Asc;
  453.     else                           /* Or -2 put object on Midheaven switch? */
  454.       j = planet[-us.objOnAsc]-MC;
  455.     for (i = 1; i <= cSign; i++)   /* If so, rotate the houses accordingly. */
  456.       house[i] = Mod(house[i]+j);
  457.   }
  458.  
  459.   /* Check to see if we are -F forcing any objects to be particular values. */
  460.  
  461.   for (i = 1; i <= cObj; i++)
  462.     if (force[i] != 0.0) {
  463.       planet[i] = force[i]-rDegMax;
  464.       planetalt[i] = ret[i] = 0.0;
  465.     }
  466.  
  467.   ComputeInHouses();        /* Figure out what house everything falls in. */
  468.  
  469.   /* If -f domal chart switch in effect, switch planet and house positions. */
  470.  
  471.   if (us.fFlip) {
  472.     for (i = 1; i <= cObj; i++) {
  473.       k = inhouse[i];
  474.       inhouse[i] = SFromZ(planet[i]);
  475.       planet[i] = ZFromS(k)+MinDistance(house[k], planet[i]) /
  476.         MinDistance(house[k], house[Mod12(k+1)])*30.0;
  477.     }
  478.     for (i = 1; i <= cSign; i++) {
  479.       k = HousePlaceIn(ZFromS(i));
  480.       housetemp[i] = ZFromS(k)+MinDistance(house[k], ZFromS(i)) /
  481.         MinDistance(house[k], house[Mod12(k+1)])*30.0;
  482.     }
  483.     for (i = 1; i <= cSign; i++)
  484.       house[i] = housetemp[i];
  485.   }
  486.  
  487.   /* If -3 decan chart switch in effect, edit planet positions accordingly. */
  488.  
  489.   if (us.fDecan) {
  490.     for (i = 1; i <= cObj; i++)
  491.       planet[i] = Decan(planet[i]);
  492.     ComputeInHouses();
  493.   }
  494.  
  495.   ciCore = ci;
  496.   return T;
  497. }
  498.  
  499.  
  500. /*
  501. ******************************************************************************
  502. ** Aspect Calculations.
  503. ******************************************************************************
  504. */
  505.  
  506. /* Set up the aspect/midpoint grid. Allocate memory for this array, if not */
  507. /* already done. Allocation is only done once, first time this is called.  */
  508.  
  509. bool FEnsureGrid()
  510. {
  511.   if (grid != NULL)
  512.     return fTrue;
  513.   grid = (GridInfo FPTR *)PAllocate(sizeof(GridInfo), fFalse, "grid");
  514.   return grid != NULL;
  515. }
  516.  
  517.  
  518. /* Indicate whether some aspect between two objects should be shown. */
  519.  
  520. bool FAcceptAspect(obj1, asp, obj2)
  521. int obj1, asp, obj2;
  522. {
  523.   int fSupp;
  524.  
  525.   if (ignorea(asp))    /* If the aspect restricted, reject immediately. */
  526.     return fFalse;
  527.   if (us.fSmartCusp) {
  528.  
  529.     /* Allow only conjunctions to the minor house cusps. */
  530.  
  531.     if ((FMinor(obj1) || FMinor(obj2)) && asp > aCon)
  532.       return fFalse;
  533.  
  534.     /* Prevent any simultaneous aspects to opposing angle cusps,     */
  535.     /* e.g. if conjunct one, don't be opposite the other; if trine   */
  536.     /* one, don't sextile the other; don't square both at once, etc. */
  537.  
  538.     fSupp = (asp == aOpp || asp == aSex || asp == aSSx || asp == aSes);
  539.     if ((FAngle(obj1) || FAngle(obj2)) &&
  540.       (fSupp || (asp == aSqu &&
  541.       (obj1 == oDes || obj2 == oDes || obj1 == oNad || obj2 == oNad))))
  542.       return fFalse;
  543.  
  544.     /* Prevent any simultaneous aspects to the North and South Node. */
  545.  
  546.     if (fSouthNode) {
  547.       if (((obj1 == oNod || obj2 == oNod) && fSupp) ||
  548.         ((obj1 == oSou || obj2 == oSou) && (fSupp || asp == aSqu)))
  549.         return fFalse;
  550.     }
  551.   }
  552.   return fTrue;
  553. }
  554.  
  555.  
  556. /* This is a subprocedure of FCreateGrid() and FCreateGridRelation().   */
  557. /* Given two planets, determine what aspect, if any, is present between */
  558. /* them, and save the aspect name and orb in the specified grid cell.   */
  559.  
  560. void GetAspect(planet1, planet2, ret1, ret2, i, j)
  561. real *planet1, *planet2, *ret1, *ret2;
  562. int i, j;
  563. {
  564.   int k;
  565.   real l, m;
  566.  
  567.   grid->v[i][j] = grid->n[i][j] = 0;
  568.   l = MinDistance(planet2[i], planet1[j]);
  569.   for (k = us.nAsp; k >= 1; k--) {
  570.     if (!FAcceptAspect(i, k, j))
  571.       continue;
  572.     m = l-rAspAngle[k];
  573.     if (RAbs(m) < GetOrb(i, j, k)) {
  574.       grid->n[i][j] = k;
  575.  
  576.       /* If -ga switch in effect, then change the sign of the orb to    */
  577.       /* correspond to whether the aspect is applying or separating.    */
  578.       /* To do this, we check the velocity vectors to see if the        */
  579.       /* planets are moving toward, away, or are overtaking each other. */
  580.  
  581.       if (us.fAppSep)
  582.         m = RSgn2(ret1[j]-ret2[i])*
  583.           RSgn2(MinDifference(planet2[i], planet1[j]))*RSgn2(m)*RAbs(m);
  584.       grid->v[i][j] = (int)(m*60.0);
  585.     }
  586.   }
  587. }
  588.  
  589.  
  590. /* Very similar to GetAspect(), this determines if there is a parallel or */
  591. /* contraparallel aspect between the given two planets, and stores the    */
  592. /* result as above. The settings and orbs for conjunction are used for    */
  593. /* parallel and those for opposition are used for contraparallel.         */
  594.  
  595. void GetParallel(planet1, planet2, planetalt1, planetalt2, i, j)
  596. real *planet1, *planet2, *planetalt1, *planetalt2;
  597. int i, j;
  598. {
  599.   int k;
  600.   real l, alt1, alt2;
  601.  
  602.   l = RFromD(planet1[j]); alt1 = RFromD(planetalt1[j]);
  603.   EclToEqu(&l, &alt1); alt1 = DFromR(alt1);
  604.   l = RFromD(planet2[i]); alt2 = RFromD(planetalt2[i]);
  605.   EclToEqu(&l, &alt2); alt2 = DFromR(alt2);
  606.   grid->v[i][j] = grid->n[i][j] = 0;
  607.   for (k = Min(us.nAsp, aOpp); k >= 1; k--) {
  608.     if (!FAcceptAspect(i, k, j))
  609.       continue;
  610.     l = RAbs(k == aCon ? alt1 - alt2 : RAbs(alt1) - RAbs(alt2));
  611.     if (l < GetOrb(i, j, k)) {
  612.       grid->n[i][j] = k;
  613.       grid->v[i][j] = (int)(l*60.0);
  614.     }
  615.   }
  616. }
  617.  
  618.  
  619. /* Fill in the aspect grid based on the aspects taking place among the */
  620. /* planets in the present chart. Also fill in the midpoint grid.       */
  621.  
  622. bool FCreateGrid(fFlip)
  623. bool fFlip;
  624. {
  625.   int i, j, k;
  626.   real l;
  627.  
  628.   if (!FEnsureGrid())
  629.     return fFalse;
  630.   for (j = 1; j <= cObj; j++) if (!ignore[j])
  631.     for (i = 1; i <= cObj; i++) if (!ignore[i])
  632.  
  633.       /* The parameter 'flip' determines what half of the grid is filled in */
  634.       /* with the aspects and what half is filled in with the midpoints.    */
  635.  
  636.       if (fFlip ? i > j : i < j) {
  637.         if (us.fParallel)
  638.           GetParallel(planet, planet, planetalt, planetalt, i, j);
  639.         else
  640.           GetAspect(planet, planet, ret, ret, i, j);
  641.       } else if (fFlip ? i < j : i > j) {
  642.         l = Mod(Midpoint(planet[i], planet[j])); k = (int)l;  /* Calculate */
  643.         grid->n[i][j] = k/30+1;                               /* midpoint. */
  644.         grid->v[i][j] = (int)((l-(real)(k/30)*30.0)*60.0);
  645.       } else {
  646.         grid->n[i][j] = SFromZ(planet[j]);
  647.         grid->v[i][j] = (int)(planet[j]-(real)(grid->n[i][j]-1)*30.0);
  648.       }
  649.   return fTrue;
  650. }
  651.  
  652.  
  653. /* This is similar to the previous function; however, this time fill in the */
  654. /* grid based on the aspects (or midpoints if 'acc' set) taking place among */
  655. /* the planets in two different charts, as in the -g -r0 combination.       */
  656.  
  657. bool FCreateGridRelation(fMidpoint)
  658. bool fMidpoint;
  659. {
  660.   int i, j, k;
  661.   real l;
  662.  
  663.   if (!FEnsureGrid())
  664.     return fFalse;
  665.   for (j = 1; j <= cObj; j++) if (!ignore[j])
  666.     for (i = 1; i <= cObj; i++) if (!ignore[i])
  667.       if (!fMidpoint) {
  668.         if (us.fParallel)
  669.           GetParallel(cp1.obj, cp2.obj, cp1.alt, cp2.alt, i, j);
  670.         else
  671.           GetAspect(cp1.obj, cp2.obj, cp1.dir, cp2.dir, i, j);
  672.       } else {
  673.         l = Mod(Midpoint(cp2.obj[i], cp1.obj[j])); k = (int)l; /* Calculate */
  674.         grid->n[i][j] = k/30+1;                                /* midpoint. */
  675.         grid->v[i][j] = (int)((l-(real)(k/30)*30.0)*60.0);
  676.       }
  677.   return fTrue;
  678. }
  679.  
  680.  
  681. /* Fill out tables based on the number of unrestricted planets in signs by  */
  682. /* element, signs by mode, as well as other values such as the number of    */
  683. /* objects in yang vs. yin signs, in various house hemispheres (north/south */
  684. /* and east/west), and the number in first six signs vs. second six signs.  */
  685. /* This is used by the -v chart listing and the sidebar in graphics charts. */
  686.  
  687. void CreateElemTable(pet)
  688. ET *pet;
  689. {
  690.   int i, s;
  691.  
  692.   ClearB((lpbyte)pet, (int)sizeof(ET));
  693.   for (i = 1; i <= cObj; i++) if (!ignore[i]) {
  694.     pet->coSum++;
  695.     s = SFromZ(planet[i]);
  696.     pet->coSign[s-1]++;
  697.     pet->coElemMode[(s-1)&3][(s-1)%3]++;
  698.     pet->coElem[(s-1)&3]++; pet->coMode[(s-1)%3]++;
  699.     pet->coYang += (s & 1);
  700.     pet->coLearn += (s < sLib);
  701.     if (!FCusp(i)) {
  702.       pet->coHemi++;
  703.       s = inhouse[i];
  704.       pet->coHouse[s-1]++;
  705.       pet->coModeH[(s-1)%3]++;
  706.       pet->coMC += (s >= sLib);
  707.       pet->coAsc += (s < sCan || s >= sCap);
  708.     }
  709.   }
  710.   pet->coYin   = pet->coSum  - pet->coYang;
  711.   pet->coShare = pet->coSum  - pet->coLearn;
  712.   pet->coDes   = pet->coHemi - pet->coAsc;
  713.   pet->coIC    = pet->coHemi - pet->coMC;
  714. }
  715.  
  716. /* calc.c */
  717.